
#ifndef BVP_H
#define BVP_H

#include <iostream>
#include <vector>
#include <stdexcept>
#include "BoundaryCondition.h"

class BVP2d
{
public:
    struct Disk
    {
        Point2d center;
        double radius;
    };

    BVP2d();

    BVP2d(Point2d center, double radius, const std::function<double(Point2d)> &f, BoundaryCondition BC);

    ~BVP2d();

    /**
     * @brief discretize the domain into grid points
     *
     * @param n number of intervals along each axis
     */
    void Discretization(int n);

    /**
     * @brief solve BVP using FD method
     *
     * @param n grid parameter n : the number of intervals along each axis
     * @return std::vector<double>
     */
    std::vector<double> Solve(int n);

    /**
     * @brief Get the Grid object
     *
     * @return std::vector<Point2d>&
     */
    std::vector<Point2d> &getGrid();

    /**
     * @brief Get the Boundary Condition
     *
     * @return BoundaryCondition&
     */
    BoundaryCondition &getBC() { return BC; }

    /**
     * @brief Get the Disk
     *
     * @return const Disk&
     */
    const Disk &getD() const { return D; }

    /**
     * @brief assign function f for solving $-\Delta u = f$;
     *
     * @param f
     */
    void assignF(const std::function<double(Point2d)> &f)
    {
        this->f = f;
    }

    /**
     * @brief print the adjacancy of each grid points
     *
     */
    void PrintAdajacancy() const;

private:
    /* data */
    // if D.radius == 0, default domian (0,1)^2
    // if not , domain (0,1)^2/D
    Disk D;

    //-\Delta u = f;
    std::function<double(Point2d)> f;

    BoundaryCondition BC;

    // store the coordinates of grid points
    std::vector<Point2d> grid;

    // store the indices of the adjacant points of one single grid points
    // Usage : adjacancy[i] = {a,b,c,d} means grid[i]'s adjacant points are grid[a],grid[b],grid[c],grid[d]
    std::vector<std::vector<int>> adjacancy;

    /**
     * @brief
     *
     * @param p
     * @return int -1 for outside points; 0 for internal points;
     *          1 for points on square boundary; 2 for points on the circle
     */
    int PointStatus(Point2d p) const;

    /**
     * @brief Internal method to link an internal point with a neiboring point
     *         if the neighboring point is outside the domain, the boundary point
     *         along the line connecting the original two points will be linked
     *         to the internal point
     *
     * @param i the index of point whose adjacancy to be updated i.e. grid[i]
     * @param h the increment along the axis of the neighboring point
     * @param axis 0 for x-axis, 1 for y-axis
     */
    inline void LinkSingle(int i, double h, int axis);

    /**
     * @brief check if the grid of internal points are connected
     *
     * @param n grid parameter n
     * @param inter_num the number of internal points in the grid
     * @return true  The Grid is Connected
     * @return false The Grid isn't Connected
     */
    bool isConnected(int n, int inter_num) const;
};

#endif